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Abstract 

Boolean satisfiability HI (/c-SAT) is one of the most studied optimization problems, as 
an efficient (that is, polynomial-time) solution to /c-SAT (for k > 3) implies efficient solu- 
tions to a large number of hard optimization problems |[2l|3l. Here we propose a mapping 
of /c-SAT into a deterministic continuous-time dynamical system with a unique correspon- 
dence between its attractors and the /c-SAT solution clusters. We show that beyond a con- 
straint density threshold, the analog trajectories become transiently chaotic iSHSHSH?], and 
the boundaries between the basins of attraction [SJ of the solution clusters become fractal 
El 13 01 5 signaling the appearance of optimization hardness ifTOl . Analytical arguments and 
simulations indicate that the system always finds solutions for satisfiable formulae even in the 
frozen regimes of random 3-SAT |11 1 and of locked occupation problems |12| (considered 
among the hardest algorithmic benchmarks); a property partly due to the systems hyperbolic 
||4j[T3 | character. The system finds solutions in polynomial continuous-time, however, at the 
expense of exponential fluctuations in its energy function. 



*The article appeared in Nature Physics 7, 966 (201 1) 
^E-mail: ercsey.ravasz @phys.ubbcluj .ro 
* E-mail: toro@nd.edu 
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Boolean satisfiability lUl (/c-SAT, A: > 3) is the quintessential constraint satisfaction problem, 
lying at the basis of many decision, scheduling, error-correction and bio-computational applica- 
tions. fc-SAT is in NP, that is its solutions are efficiently (polynomial time) checkable, but no 
efficient (polynomial time) algorithms are known to compute those solutions 0. If such algo- 
rithms would be found for fc-SAT, all NP problems would be efficiently computable, since fc-SAT 
is NP-complete EII3. 

In k-SAT there are given N Boolean variables {xi, . . . , x^}, xi G {0, 1} and M clauses 
(constraints), each clause being the disjunction (OR, denoted as V) of k variables or their negation 
(x). One has to find an assignment of the variables such that all clauses (called collectively as 
a formula) are satisfied (TRUE =1). When the number of constraints is small, it is easy to find 
solutions, while for too many constraints it is easy to decide that the formula is unsatisfiable 
(UNSAT). Deciding satisfiability, in the 'intermediate range', however, can be very hard: the 
worst-case complexity of all known algorithms for fc-SAT is exponential in N. 

Inspired by the mechanisms of information processing in biological systems, analog comput- 
ing received increasing interest from both theoretical [14l[l5l[i6| and engineering communities 
|[T7l[T8l[T9l[20l[2Tl. Although the theoretical possibility of efficient computation via chaotic dy- 
namical systems has been shown previously 1 15 1, nonlinear dynamical systems theory has not been 
exploited for NP-complete problems in spite of the fact that, as shown by Gu et al.| 19 |, Nagamatu 
et al. h20il and Wah et al. [21], fc-SAT can be formulated as a continuous global optimization 
problem |[T9l , and even cast as an analog dynamical system ll20l[2T1l . 

Here we present a novel continuous-time dynamical system for fc-SAT, with a dynamics that 
is rather different from previous approaches. Let us introduce the continuous variables I19II si G 
[—1,1] , such that 5^ = — 1 if the z-th variable {xi) is FALSE and 5^ = 1 if it is TRUE. We define 
Cjni = 1 for the direct form (x^), Cmi = —1 for the negated form (x^), and Cmi = for the absence 
of the z-th variable from clause m. Defining the constraint function Km{s) = Ililill ~ 
CmiSi) corresponding to clause m, we have Km G [0, 1] and Km = if and only if clause m is 
satisfied. The goal would be to find a solution s* with 5* G { — 1, 1} to £^(s*) = 0, where E 
is the energy function E{s) = J2m=i ^m{s)'^' If such s* exists, it will be a global minimum 
for E and a solution to the fc-SAT problem. However, finding s* by a direct minimization of 
E{s) will typically fail due to non-solution attractors trapping the search dynamics. In order 
to avoid such traps, here we define a modified energy function V{s^a) = Ylm=i ^mKmi^Y^ 
using auxiliary variables am ^ (0, oo) similar to Lagrange multipliers |[20l|2Tl. Let us denote by 
%N the continuous domain [—1, 1]^. Its boundary is the A^-hypercube Qn = OT-Ln with vertex 
set Vn = {— l,!}^ C Qn' The set of solutions for a given fc-SAT formula, called solution 
space, occupies a subset of Vn- Solution clusters are formed by solutions that can be connected 
via single-variable flips, always staying within satisfying assignments [22J. Clearly, V > in 
= T-Ln X (0, oc)^, and V{s^ a) = within Vn if and only if s = s* G Vn is a fc-SAT solution, 
for any a G (0, oc)^. We now introduce a continuous-time dynamical system on Q through: 



dsj 
dt 



M 



{-VsV{s,a)). ^ 2amCmi 



Kmi{s)Km{s) , Z = 1, . . . , 



(1) 



dar 



"m 



amKm{s) 



m = 1, . . . ,M 



(2) 



dt 
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where V5 is the gradient operator with respect to s, and Kj^i = i^m/(l — CjniSi). The initial 
conditions for s are arbitrary s(0) ^T-Ln, however, for a they have to be strictly positive, a^(0) > 
(e.g., am{0) = 1). The fc-SAT solutions s* G Vn are fixed points of ([T]|2]),/(9r any a E (0, oc)^. 
The fc-SAT solution clusters are spanning piecewise compact, connected sets in Qat, and every 
point in them is a fixed point of ([T]|2]) (Supplementary sect. A). System ([T]|2]) has a number of 
key properties (see Supplementary Information), (i) The dynamics in s stays confined to T-L^. (ii) 
The fc-SAT solutions s* G Vn are attractive fixed points of ([l]|2]). In particular, every point s 
from the orthant of a fc-SAT solution s* with the property |sp > iV - 1 + (/c - + 1)^ is 

guaranteed to flow into the attractor corresponding to s*. (iii) There are no limit cycles, (iv) For 
satisfiable formulae the only fixed point attractors of the dynamics are the global minima of V with 
y = 0. Note that in principle, the projection of the dynamics onto T-Ln could be stuck in some 
point s, while da/dt ^ indefinitely. This does not happen here, as shown in Supplementary 
sect. E. Moreover, analytical arguments supported by simulations indicate that the trajectory will 
leave any domain that does not contain solutions, see the discussion in Supplementary sect. El. 
Note, that the constraint functions (hence their satisfiability) depend directly only on the location 
of the trajectory in Km = Km{s), and not on the auxiliary variables. The dynamics in the 
a-space is simple expansion, and for this reason the features of the full phase space lie within its 
projection onto T^tv- One can actually eliminate entirely the auxiliary variables from the equations 
by first solving (2) to give a^(t) = a^(0) exp ^ Km{s{T))dT^ then inserting it into ( 1 ). 

Another funaamental feature of ([T][2]) is that it is deterministic: for a given formula /, any 
initial condition generates a unique trajectory, and any set from T-L^ has a unique preimage arbi- 
trarily back in time. Hence, the characteristics of the solution space are reflected in the properties 
of the invariant sets ||7| of the dynamics ([T][2]) within the hypercube T-Ln. The deterministic na- 
ture of ([T]|2]) allows us to define basins of attractions of solution clusters by colouring every point 
iuT-Ln according to which cluster the trajectory flows to, if started from there. These basins fill 
AT up to a set of zero (Lebesgue) measure, which forms the basin boundary 0, from where the 
dynamics (by definition) cannot flow to any of the attractors. A fc-SAT formula / can be repre- 
sented as a hypergraph Q{f) (or equivalently, a factor graph) in which nodes are variables and 
hyperedges are clauses connecting the nodes/variables in the clause. Pure literals are those that 
participate in one or more clauses but always in the same form (direct or negated); hence they 
can always be chosen such as to satisfy those clauses. The core of G{f) is the subgraph left after 
sequentially removing all the hyperedges having pure literals [23|. For simple formulae (such as 
those without a core), the dynamics of ([T]|2]) is laminar flow and the basin boundaries form smooth, 
non-fractal sets (Fig|T^,c and Figj2]top two rows). Adding more constraints G{f) develops a core, 
the spin equations ([T]) become mutually coupled, and the trajectories may become chaotic (Fig|Tj3, 
Supplementary sect. F, Fig. 12) and the basin boundaries fractal |[7l[8l[9l (Fig. [TJl, Fig. [2} Fig. 8). 
Therefore, as the constraint density a = M/N is increased within predefined ensembles of formu- 
lae (random fc-SAT, occupation problems, fc-XORSAT, etc.) a sharp change to chaotic behaviour 
is expected at a chaotic transition point a^, where a chaotic core appears with non-zero statistical 
weight in the ensemble as ^ oc. As an example, let us consider 3-XORSAT. In this case, due 
to its inherently linear nature, it is actually better to work directly with the parity check equations 
as constraints, instead of their CNF form. The chaotic core here is a small finite hypergraph, and 
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Figure 1: Chaotic behaviour. Five, closely started sample trajectories projected onto (si, 52, 53) 
a), for a 3-SAT formula with N = 200, a = 3 and b) for a hard formula, N = 200, a = 
4.25. The colour indicates the energy E (colour bar) in a given point of the trajectory. While 
for easy formulae the trajectories exhibit laminar flow, for hard formulae they quickly become 
separated, showing a chaotic evolution. Taking a small 3-XORSAT instance with = 15 (see 
Supplementary sect. G) we fix a random initial condition for all Si, except si and 52 which are 
varied on a 400 x 400 grid and we colour each point according to the solution they flow to for c) 
7 = 0.6 (instance shown on Fig. 8e) and d) for 7 = 0.8 (instance shown in Fig. 8f)). 
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thus OL^ coincides with the so-called dynamical transition point computed exactly by Mezard et 
al. |[24l (see Supplementary sect. G and Fig. 8). Note, a core can be non-chaotic, and thus the 
existence of a core is only a necessary condition for the appearance of chaos and in general the 
two transitions might not coincide. Further increasing the number of constraints (within any for- 
mula ensemble) unsatisfiability appears at the threshold value ols > beyond which almost all 
formulae are unsatisfiable (UNSAT regime)!!!! [!2 IM ESI EZl Ell. The closer a is to 
as, the harder it is to find solutions, and beyond the so-called freezing transition point af < ag 
(called the frozen regime) all known algorithms take exponentially long times or simply fail to 
find solutions ||!!ll!3. A variable is frozen if it takes on the same value for all solutions within 
a cluster, and a cluster is frozen if an extensive number of its variables are frozen. In the frozen 
regime all clusters are frozen and they are also far apart ( 0{N) Hanmiing distance) |[!!l[!2l. For 
random 3-SAT (clauses chosen uniformly at random for fixed a) ag = 4.26 L27]| . o^/ = 4.25 |[28ll 
and all known local search algorithms become exponential or fail beyond a = 4.21 |[29l , while 
Survey Propagation ll25l based algorithms fail beyond a = 4.25 |[28ll . Since the frozen regime 
is very thin in random 3-SAT, Zdeborova and Mezard [12] have introduced the so-called locked 
occupation problems (LOPs). In LOPs all clusters are formed by exactly one solution, hence they 
are completely frozen and the frozen regime extends from the clustering (dynamical) transition 
point Id to the satisfiability threshold Ig, and thus it is very wide llT2ll . An example LOP is random 
+l-in-3-SAT [12], made of constraints that have no negated variables and a constraint is satis- 
fied only if exactly one of its variables is 1 (TRUE). In +l-in-3-SAT Id = 2.256, Ig = 2.368, 
and beyond Id all known algorithms have exponential search times or fail to find solutions (here 
/ = 3M/N). 

As chaos is present for satisfiable formulae, that is, when system ([T]|2]) has attracting fixed 
points, it is necessarily of transient type. Transient chaos S [51 |6l [3 is ubiquitous in systems 
with many degrees of freedom such as fluid turbulence |[30l . It appears as the result of homo- 
clinic/heteroclinic intersections of the invariant manifolds of hyperbolic (unstable) fixed points of 
([T]|2]) lying within the basin boundary Il7l[8l[9l, leading to complex (fractal) foliations of the phase 
space (see Supplementary sect. F). We observed the prevalence of transient chaos in the whole 
region < a < ag for all the problem classes we studied. Interestingly, the velocity fluctuations 
of trajectories in the chaotic regime are qualitatively similar to those of fluid parcels in turbulent 
flows as shown in Supplementary sect. K. Our findings suggest that chaotic behaviour may be a 
generic feature of algorithms searching for solutions in hard optimization problems, corroborating 
the observations by Elser et al. lfTOl using a heuristic algorithm based on iterated maps. 

In the following we show results on random 3-SAT and +l-in-3-SAT formulae in the frozen 
regime, however, the same conclusions hold for other ensembles that we tested. To investigate 
the complexity of computation by the flow ([l]|2]), we monitored the fraction of problems p{t) not 
solved by continuous time t, as function of N and a. Figs. |3^,c show that even in the frozen 
phase, the fraction of unsolved problems by time t decays exponentially with t, that is, by a law 
p{t) = j,^-KN)t^ jj^g decay rate A(A^) obeys A(A^) = bN~^, with /3 2^ 1.6 in both cases, see Fig. 
[3|3,d. From these two equations, the continuous-time N) needed to solve a fixed (1 — p)ih 
fraction of random formulae (or to miss solving p-th fraction of them) is: 

t{p,N) = b-^N^ln{r/p) (3) 
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Figure 2: Attractor basins. For a random 3-SAT instance with = 50 we vary a by successively 
adding new constraints. Fixing a random initial condition for Si, z > 3, we vary only si and 52 
on a 400 x 400 grid, and we colour each point according to the solution (first column) or solution 
cluster (second column) they flow to. Each colour in a given column represents a solution or 
solution cluster respectively, however colours between columns are independent. The third column 
represents the analog search time t needed to find a solution (see colour bar) starting from the 
corresponding grid point. Maps are presented for values of a = 3.5, 3.7, 3.9, 4.1, 4.16, 4.2, 4.24. 
Easy formulae are characterized by smooth basin boundaries and small search times. Note that we 
only see the solutions (and clusters) that reveal themselves in the (51, 52) plane, others might not 
be seen. For hard formulae the boundaries and the search time maps become fractal. 
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Figure 3: Computational complexity properties, a) The fraction of problems p{t) not yet 
solved by continuous-time t for 3-SAT at a = 4.25, for N = 20, 30, 40, 50, 60, 80, 100, 125, 150 
(colours). Averages were done over 10^ instances for each N. For each instance the dy- 
namics was started from one random initial condition. Black continuous lines show the de- 
cay p{t) = r exp(— A(A^)t). b) The decay rate follows A(A^) = bN~^, with /3 :^ 1.66. 
c) The fraction of problems p{t) unsolved by time t for +l-in-3-SAT at / = 2.34, for N = 
20, 25, 30, 35, 40, 50, 60, 70, 80. For each instance the dynamics was started in parallel from 10 
random initial conditions, averages were taken over 10"^ instances for each N. Black continuous 
lines show the same exponential decay as in a), d) The decay rate shows the same behaviour as in 
b) with exponent: /3 1.68. 
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indicating that the continuous time (CT) needed to find solutions scales as a power-law with N. 
Eq. ^ also implies power-law scaling for almost all hard instances in the ^ oc limit (Sup- 
plementary sect. H). The length in 7^ at of the corresponding continuous trajectories also scales 
as a power-law with N (Supplementary Fig. lib, sect. J). However, note that this does not 
mean that the algorithm itself is a polynomial-cost algorithm, as the energy function V can have 
exponentially large fluctuations. As the numerical integration happens on a digital machine, it 
approximates the continuous trajectory with discrete points. Monitoring the fraction of formulae 
left unsolved as function of the number of discretization steps Ugtep in the frozen phase, we find 
exponential behaviour for rigtepip^ N) (Supplementary sects. I,J, Fig. 10). The difference between 
the continuous- and discrete-time complexities is due to the wildly fluctuating nature of the chaotic 
trajectories (see Fig. [TJ3 and Methods) in the frozen phase. Compounding this, we also observe the 
appearance of the Wada property [TtIISI in the basin boundaries. Fig. |4j A fractal basin boundary 
has Wada property if its points are simultaneously on the boundary of at least three colours/basins. 
(An amusing method that creates such sets uses four Christmas ball ornaments [7].) Although the 
Wada property does not affect the true/mathematical analog trajectories, owing to numerical er- 
rors, it may switch the numerical trajectories between the basins. Since the clusters are far (0{N)) 
apart, the switched trajectory will flow towards another cluster into a practically opposing region of 
T-Ln until it may come close again to the basin boundary etc., partially randomizing the trajectory 
in T-Ln- 




Figure 4: Wada property. Basin boundaries are shown for +l-in-3-SAT for an instance at = 
30 and / = 2.28. Fixing a random initial condition for all Si, i > 3 we vary only si and 52 on a 
200 X 400 grid and colour the points according to three different solutions they flow to. Successive 
magnifications illustrate the Wada property: the points on the basin boundaries are simultaneously 
on the boundary of all three basins implying that large enough magnifications will contain all three 
colours (although the blue-green boundary seems void of red in the third panel, panels four and 
five show that red is actually present). 

We conjecture that the power-law scaling of the continuous search times ([3]) is due in part to 
a generic property of the dynamical system ([T]|2]), namely that it is hyperbolic ||4l[6l[T3l, or near- 
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hyperbolic. It has been shown that for hyperbolic systems the trajectories escape from regions far 
away from the attractors to the attractors at an exponential rate, for almost all initial conditions 
im [6l [131. That is, the fraction of trajectories still searching for a solution after time t decays 
as e~^^(Supplementary Fig. 13), where n is the escape rate. Thus, /^~^ can be considered as a 
measure of hardness for a given formula. When taken over an ensemble at a given a, this property 
generates the exponential decay for p{t) with an average escape rate A. 

The form of the energy function V incorporates the influence of all the clauses at all times, and 
in this sense the system ([T][2]) is a non-local search algorithm. As shown before, the auxiliary vari- 
ables can be eliminated, however, they give a convenient interpretation of the dynamics. Namely, 
one can think of them as providing extra dimensions along which the trajectories escape from local 
wells, and their form ([2]) provides positive feedback that guarantees their escape. Clearly, these 
equations are not unique, and other forms based on the same principles may work just as well. 

Methods 

To simulate ([l]|2]), we use a 5-th order adaptive Cash-Karp Runge-Kutta method with monitoring 
of local truncation error to ensure accuracy. In order to keep the numerical trajectory within a 
tube of small, preset thickness around the true analog trajectory in (Supplementary Fig. 9), the 
RK algorithm occasionally performs an exponentially large number of discretization steps nstep- 
However, this only happens for hard formulae, when the analog trajectory has wild, chaotic fluctu- 
ations. For easy formulae both p{t) and p{nstep) decay exponentially as shown in Supplementary 
Fig. 10a, inset. 
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A Supplementary Information 

Here we provide derivations and discussions about the properties of the dynamical system ([T][2]) 
and additional supporting figures and text. 

Recall from the main text the definitions: 

M 

y(a,s) = 5^a^K^, (4) 

m=l 

where 

N 

Km = 2-^l[{l-CmjSj) . (5) 
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Note that in fc-SAT there are at most k terms in the product above, hence < < 1 for all m. 
The system of ODEs (fTpl) defined in the main document is: 



5i = ^ = --^V{a, 5) = ^ 2amCmiKmiKm , z = 1, . . . , TV , (6) 

am = = CimKm , m = 1, . . . , M , (7) 

at 

with a^(0) > 0, m = 1, . . . , M ( for example a^(0) 1), and where 

i^m. = 2-^ f] (1 - C^,5,) = _ . (8) 



A. Free variables, solution clusters and attractors. Solution clusters are defined by solutions 
that can be connected via single-variable flips, always staying within satisfying assignments. For 
example, consider two k-SAT solutions that differ in exactly one variable, let's say in sj (at Ham- 
ming distance of 1), thus forming a solution cluster of two points in Vn = { — ^A}^- Then any 
point on the sj axis in the [—1,1] continuous segment is a fixed-point (s = 0, d = 0) of the dy- 
namics (6]|7 ), because for the two solutions the value of sj is irrelevant (called a "free variable"), 
all the clauses being satisfied by the other variables taking values ±1. For our dynamical system, 
the corresponding attractor is an edge of the A^-hypercube Qat, a compact, connected domain/set. 
It is certainly possible that several solutions in Vn be connected via single variable flips, with two 
examples given in Figs. [5^-b. Thus, the attractors of the dynamical system ([6]|7]) (or (1) in the 
main text) in general are compact, connected sets of Qn spanned by the k-SAT solution clusters 
and by only those (see also sections D and E). Every point from this attractor set is a fixed-point 
of the continuous-time (CT) dynamics. 



B. Spin variables s remain in 7/ at = [-1, 1]^. From Eq. 

M M 



Si = 2 ^ amCmi{l - CmiSi)K^i = ^ ^ CLfuK^i [(1 " 5i)5c^,,l " (1 + Si)6c^^-i] (9) 



m=l m=l 



This dynamics keeps spins within [—1,1], because at s^ = 1 (si = —1) we have Si < (si > 0) 
for any Cmi ^ {-1, 1}. 



C. Stability of all fc-SAT solutions. Domains of attraction. As seen in the main text, all the 
s points for which y = 0, are fixed points of the dynamics ([6]|7]). As shown in Section A, some 
of these points are not necessarily from Vn, they can be from the compact domain corresponding 
to a cluster of solutions from Vn, when there are free variables. We prove stability by showing 
that in a vicinity of non-zero volume in at of a fc-SAT solution, the function R = J2iLi ^1 
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Figure 5: Example attractors spanned by solution clusters. Two examples with = 3, = 3, 
M = 3, see the figure for the sets of clauses in each case. In a the solution cluster is the set 
of vertices {A, O, C, D}, while for the dynamical system the corresponding attractor is the 
continuous union of segments AO U OB U BC U CD (and the equivalent in the s-space). In b, 
the solution cluster is formed by the four vertices of the shaded square-plate tt and vertex A, while 
for the dynamical system the corresponding attractor is tt U AO. 



monotonically increasing until the dynamics reaches a point in s-space for which V = (hence 
all Km = 0), that is, the trajectory reaches the cluster's domain. From ([6]): 

N M N 

R = 2 ^ ^ SiSi — 4 ^ ^ djYiKjYi ^ ^ CjYiiSiKjYii . (10) 

i=l m=l i=l 

Assume that s* G Vn is a fc-SAT solution. Recall that we have k variables in each clause (all 
the rest have Cmi = 0). Choosing a clause Km, let us denote the indices of the variables included 
as z = 1, . . . , A:, and order them such that the variables j = 1, . . . ,p satisfy the clause, meaning 
that (1 — CmjSj) = 0, and the variables /=p+l,...,A:do not satisfy the clause, and thus 
(1 — CmiSi) = 2. Clearly, such a 1 < p < A: always exists, since = 0. If we are in the 

corresponding A^-dimensional octant (orthant) of s* then sgn{si} = 5* for all s from this orthant 
with < l^^l < 1. Then CmjSj = 1 implies CmjSj = \sj\ for all j = 1, . . . ,p and = — 1 

implies Cmisi = — for all / = ]9 + 1, . . . , A: in this orthant. Thus: 

Zm = ^ CmiSiKmi ^ Km^ _™ ' ^ ^rn _\ \ - ^ , , (H) 

i=i i=i ^ \j=i ^ i^^i i=p^i ^ + i^^iy 

Clearly, for any a > with a < miuj \sj\, we have — \sj\) > cr/(l — a) (for all j) and 

+ < 1/2 (since < 1). Therefore 

Z.>K.(^-^). (12) 
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Hence, if a > ^ , then Zm > (for all m) and from ( 10) i? = J2m ^^mKmZm > 0, that is 

s* = (l,l,s*,...) 




k-l 
' k + l 



Figure 6: Domain of attraction, 2D illustration. If s* G VV is a k-SAT solution, and sgn{si) = 

5*, > (A: — + 1) are satisfied for alH = 1, . . . , (inside the blue square) then dR/dt > 
0. The dynamics cannot leave from inside the domain CAr(s*) (marked with red) as \si\ < 1 
(see section B). It follows that R grows until it hits the domain of the solution cluster attractor 
containing s*, or s* itself, meaning that CAr(s*) is part of the attractor's basin. 

the trajectory is strictly increasing its distance from the origin (unless all = 0, which means 
the trajectory is on the attractor, where R = 0). Since {k — p)/{k + p) is a decreasing function of 
p, we can choose p = 1 and set 

k-l 

o' = 

fc + l 

to define a vicinity of s*, with l^^l > a, i = 1, . . . , within which we are guaranteed R > 0. 
Next, consider the corner domain Cat (s* ) around s"^ of Hn cut away by the A^-dimensional sphere 
of radius \/A^ — 1 + centered in the origin. See Fig. |6|for an illustration in 2D. Trajectories in 
all points within this domain have the property i? > or i? = 0, in the latter case the point being 
on the attractor. For points with i? > 0, the trajectory must necessarily flow towards the boundary 
part {Qn) of the hypercube of this domain CAr(s*), away from the surface of the sphere until it hits 
s*, or the attractor of the solution cluster of which s* is part of (in case of free variables), lying 
within CAr(s*). Clearly R cannot increase beyond A^, in which case the trajectory is at s*. 



D. There are no limit cycles in s. Having a limit cycle in s which is not a fixed-point means 
that Si{t) is a periodic function of t of period T > 0, and consequently all its derivates, including 

s,{t) = s,{t + nT) , s,{t) = s,{t + nT) (13) 

for alH = 1, . . . , A^ and all integers n. Since the Kj^ are functions of time only through the si 
variables, this implies that all Kj^ are periodic functions of time as well. From ([6]) it follows: 

M 

S^{t + nT)=Yl ""rnit + nT)2CrmKm{s{t))Krm{s{t)). (14) 

m=l 
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Formally, the solution for from ^ can be written as 



cim{t) (to) exp 



f drKMr)) 



(15) 



Using the periodicity of Km this leads to a^(t + nT) = ajn{t)e^^'^ with = Km{s{r))dr. 
From ([14]) and ([13]) it follows: 

M 



J2 ^mit) (e^^-^ - 1) 2CrmKm{t)Krm{t) = , Z = 1, . . . , TV 



m=l 



which has to hold for any integer n and all times t. Every sum in ( |T6| ) is of the type 

5^ - 1) /mz = , V n G Z , 



(16) 



(17) 



Where Xm = e^^. Since Km G [0, 1] it follows that /m ^ and thus Xm > 1- Assume that the 
XjYi~^ are all different from each other, and consider x^* to be the largest of them. Then, clearly. 



increasing n without limit, all other terms in the sum of ( [T7| ) become arbitrarily small compared 
to the m* term and since x^* > 1, this forces fm*i = 0. Thus the m* term must be absent from 
the sum ( [T7] ), and we can repeat the procedure with the next largest x^ term leading to fmi = 0, 
etc. Hence, all fmi = and thus Si = J2m f'rni = 0. If the x^-s are equal "in blocks", with a 
similar procedure we can show that for every block h we must have {x^ — 1) X]mG6 f'rni = 0, for 
all n. This implies that either = 1, or ^m^h f'rni = 0. In the former case we thus must have 
for all m e b, Im = and hence Km = (recall, that Km > 0) which actually implies that 
the corresponding fmi = 0, since fmi oc Km- Hence, either way, J2meb f™ ^ ^ blocks, 
and therefore, again, Si = fmi = meaning that we are in a fixed point, contradicting our 
original assumption. 



E. Fixed-points and attractors in T^tv. We have seen that the dynamics in Q admits fixed point 
solutions (s = 0, d = 0) only on the boundary Qn, either as isolated points from Vn (fc-SAT 
solutions), or in form of compact, connected domains (if there are free variables) corresponding 
to fc-SAT solution clusters. This holds because d = if and only if Km = for all m, which is 
possible (due to their product form, (jSj) only if some of the spins are ±1, hence they form fc-SAT 
solutions. The question is whether there are other stable fixed points within %n (thus in s-space), 
but with not all Km = 0, in which the dynamics could get stuck indefinitely. In principle it could 
happen that there are points s G T-Ln for which s = but the dynamics ^ of the auxiliary 
variables is not able to unstuck the spin variables from there. The answer is negative that is, there 
are no such stable fixed points and here we sketch its proof. However, as we will see, there can 
be unstable fixed-points, which play an important role in the dynamics within T^at. In subsection 
E.l below we discuss how the trajectory must leave any domain V within T-Ln that does not have 
a solution in it. 

Let us assume that the dynamics at time t arrives into a point s = s(f) G T^at for which: 



kit) 



d_ 

ds, 



V{a,s] 



M 

2 5] a, 

m=l 



'miKm{^^Kmi{^^ — 



1, 



(18) 
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and not all K^-s are zero. Let us denote = am{t), = Km{s). Then, using dsl), system 



( 18 ) can be thought of as a homogeneous system of equations for the variables: 



M M 

^m^m -i _ ^ a^^^mi = , Z = 1, . . . , TV . (19) 

m=l ^ ^^^^^ m=l 

As there are more clauses than variables (M > N), system ([19]) in general defines an M — r- 
dimensional domain in a-space, namely, the left null-space ker(f7^) of the M x matrix U = 
{umi} of rank r < N. Since not all K^-s are zero, the variables will change over time, 
according to ([7]). In order for the dynamics to be stuck in s indefinitely, one must have condition 

|l9| ) hold for all later times t + r. Using a^(t) = a^e^^^^^^^^ = dme^^"^, the condition 
becomes: 

M 

e'^^cimUmi = , z = 1, . . . , TV , for all r > . (20) 

m=l 



If there is a single rrf for which > for all m, then clearly, the Ihs of ( |20| ) becomes 
dominated by this maximum exponential term and ( [20| ) becomes violated, and hence the dynamics 
becomes unstuck from s. The only way that ( [20| ) has a chance to hold for arbitrary r, if the 
K^-s are equal in blocks of at least size two (since Umi oc Cmi G {—1,0, 1}, cancellations 
are possible). That means that there are at least M/2 equations expressing the equalities of the 
corresponding Kj^-s. Since there are N spin variables and we are looking at a = M/N > 2 
(no point solving SAT in the very easy phase) there are more (nonlinear) equations than variables 
and thus it drastically reduces the chances of finding an s solution satisfying these. Assuming 
that such s exists (one cannot exclude it in general), we need to analyze this case, in particular 
whether the equality in blocks of the Km-^ is an "attractive" condition by the dynamics for some 
s. In other words, we need a stability analysis of s. It is important to note that s is not a standard 
fixed point in the traditional sense. Since the dynamics is also driven by the auxiliary variables, we 
have the continuous domain of the left null-space ker(C7^) as "fixed point" for the dynamics. As 
soon as the variables leave this domain, the dynamics gets unstuck from s, and only then. For 
this reason, the linear stability analysis around s is somewhat different from a standard stability 
analysis. 



Thus, let us assume that ( [20| ) holds for arbitrary r for some s. Since the dynamics is in the 
(s, a) -space, let us consider a small deviation e such that a = a + e, > 0, m = 1, • • • , M, 
|e| <^ 1, at T = 0. Clearly, one can easily choose the Cm variables, such that 

M 

Y^mUmiT^O^ i = l,...,A^, (21) 

m=l 

that is, e is any vector that is not from the left null-space of the {umi} matrix. For example, if 



uii ^ 0, then e = 0, . . . , 0), 5 > 0, 5 <C 1 would suffice. From (18) it then follows, 

that the Si variables start changing, as 5^ = 2 J2m ^mUmi 7^ 0. Recall that in fc-SAT we are only 
considering problems in which every variable is present in both its direct and negated forms, in 
different clauses. If this wasn't the case, that is, if a variable Si when present in a clause would be 
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only present in its direct form (no negated form anywhere), then we could easily set the value of 
that variable to +1 and automatically satisfy all the clauses in which it is present and reduce the 
problem to a smaller one. Additionally, if all c^^-s would have the same sign, one could never 
have ( [T9| ) satisfied in the first place, since > and 1 — CmiSi > 0. This means, that for 
any Si ^ 0, there will be clauses for which the Km-s will increase and others for which it will 
decrease. Since these variations are continuous, the e shift vector can always be chosen such that 
there will appear a single largest clause Kj^* , resulting in the case already discussed above, with 
a leaving exponentially fast the left null-space ker(C7^) showing that the dynamics is unstable in 
s. 

Clearly, it is in general mathematically possible for the Kj^-s to be equal in blocks such that 
( [20| ) holds for all r, defining regions of zero (Lebesgue) measure in T-L^. As we have just shown, 
these fixed regions or "points", however, are all unstable for the dynamics. The linearized neigh- 
borhoods of these points are characterized by the stable and unstable subspaces spanned by the 
eigenvectors with contracting and expanding eigenvalues, respectively Bll . For this reason, these 
type of fixed points are called hyperbolic fixed points or saddles. Moving away from the linearized 
neighborhoods, these two subspaces form the stable and unstable manifolds of the saddle, which 
then extend endlessly into the phase space. While two stable and two unstable manifolds can 
never cross, the stable and unstable manifolds can intersect, forming either homoclinic intersec- 
tions (when the two manifolds belong to the same saddle) or heteroclinic intersections (they come 
from two different saddles), see Ref lISTTl . Chaotic dynamics appears as the result of homoclinic 
or heteroclinic intersections. 

E.l Escape from an arbitrary domain. So far we have shown an important property, that is, the 
dynamics cannot be captured by a non- solution J point (for any formula). We have also shown 
in Section D that it cannot be captured by a limit cycle either. The question remains, however, 
whether it could be captured by some other type of non-solution attractor, possibly even a chaotic 
attractor. One way to show that this cannot happen, is to prove that the trajectory cannot stay 
confined for arbitrarily long times within an arbitrary domain V in the s-space (V C T-Ln) that 
does not contain a solution. The following description is a brief sketch for what happens during 
the dynamics, and it forms the elements of a proof (to be published elsewhere). 

Let us assume that at some point in time to, the trajectory is in V, s{to) E V. After eliminating 
the auxiliary variables from ([6]) using their expression from ( [TS] ) , the rate of change for the spin 
variables becomes: 

m=l 

As s is confined to V, there is always a subset of constraints (Kj^-s) that are bounded away 
from zero within V (otherwise we are in a solution, which contradicts the assumption that V 
has no solutions in it). After long enough times, the exponentials in ( |22| ) involving these con- 
straints grow very large, and unless they balance each other perfectly (due to the c^^-s) in ev- 
ery s point along the trajectory in V, one of them will overtake the others and overgrow them 
(exp J^^ drKjj^ > e^(^~^o)), at least for a while. In this time interval the sign for some of the 
dsi/dt-s will stabilize into either just positive or just negative and their magnitude will grow expo- 
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nentially fast {\Ae^'^ — Be^^\ ^ gmax{a,6}t^ When that happens, at large enough t, the trajectory 
will outburst from V (a finite domain). It can also happen, however, that the constraints oscillate 
around each other, and hence their role as to which one is the leading one in (22), changes over 
time. This can result in an oscillatory behavior in the corresponding dsi / dt-s and of the trajectory. 
However, even in this case, all terms corresponding to the constraints are growing exponentially 
fast, and the differences between the growing exponentials have increasingly wilder fluctuations, 
which eventually leads to the outburst of the trajectory from V. One can show that any situation in 
which the exponential terms corresponding to the constraints in all the N equations perfectly can- 
cel each other (so that there is no leading exponential) will be unstable against small perturbations, 
similarly to the case of a fixed point. 



F. Finite Size Lyapunov Exponents (FSLE). We have employed the FSLE method ||32l from 
nonlinear dynamics theory to provide a distributed measure of chaos in our system. The FSLE 
describes the local average strength of exponential separation of trajectories and it has been exten- 
sively used to analyze turbulent flows and processes, including in the atmosphere 1 33 ] and oceans 
lf34ll . The FSLE in a point s is given by (j){s, sq, ^) = Ine/eo), where eo = \s — is a small 
initial separation of two points with chosen in a random direction, and r is the time needed for 
the separation of the corresponding trajectories started from these two points to reach the given 
separation s. The average (•) is over the random directions of s^ In FigjT]^ = 30^o and is 
averaged for 50 different randomly oriented initial separations. 



a=3.00 b a=4.25 




Figure 7: Finite-Size Lyapunov Exponents. The colour maps show Finite Size Lyapunov Expo- 
nents measured on a problem with A: = 3, = 50 at a: o: = 3 and b: o: = 4.25. Initially, all sj 
have a random, but fixed value except si and 52 which are varied along a 400 x 400 grid. FSLE 
values are coded with colors, shown by the color bar. Enhanced chaotic behavior appears for hard 
formulae. 
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G. Random fc-XORSAT. In fc-XORSAT there are given M parity check equations (constraints) 
|[35l|36l[37l, each involving k specified Boolean variables and a parity bit = {0, 1}: 

+ . . . + = ymij^oA 2), m = 1, .., M. (23) 

The k distinct variables present in a parity check equation are chosen uniformly at random from the 
set of xi, . . . , xat variables. The goal is to assign Boolean, {0, 1} values to the variables such that 
all parity checks are satisfied. As shown in propositional calculus, any propositional formula can 
be transformed into its conjugate normal form (CNF), and hence boolean decision problems can all 
be cast into a fc-SAT problem. Once fc-XORSAT is transformed into fc-CNF (Conjugate Normal 
Form), it becomes a fc-SAT problem over a specific ensemble of clauses, hence our dynamical 
system (1) can be used to solve fc-XORSAT (or any boolean decision problem). For = 3 we use 
four 3-SAT clauses to encode one parity check equation. 

Random XORSAT also goes through phase transitions |[36l [37l when increasing the con- 
straint density 7 = M/N. The dynamical phase transition takes place at (for 3-XORSAT 
7(i = 0.8185), when the unique solution cluster existing at 7 < 7^^ breaks into an exponentially 
large number of clusters. The second phase transition is the SAT/UNSAT transition at 7c > 7(i, 
with 7c = 0.9179 for 3-XORSAT |[37]| . Constraint satisfaction formulae can be represented as hy- 
pergraphs with nodes representing the variables and with hyperlinks representing the constraints 
connecting the variables present in them. To find the various transitions for fc-XORSAT, however, 
it is actually better not to bring it into CNF form (we do that for simulations with our system 
([6]),(|7]). Using this representation, it was proven |37 1 that corresponds to the constraint density 
where the hyperloops (in this non-CNF form hypergraph) appear with non-zero statistical weight 
in the limit of ^ oc. As chaotic behavior in our dynamical system (1) appears already at 
one of the smallest (finite) hyperloop motifs, it therefore also appears at exactly the same that 
Mezard et.al. calculated |[37ll . Fig [8^) shows the full hypergraph for a small specific instance of 
3-XORSAT without a core at 7 = 9/14 = 0.64 (after performing the leaf-removal algorithm 
described in 1 37 1 the remaining hyperloops form the core of the XORSAT instance). As the graph 
in Fig [8^) has no core (no hyperloops), the corresponding basin boundary is indeed smooth, i.e., 
there is no (transient) chaos, also illustrated in Fig [8^). However, once we add more constraints 
(for example the specific ones shown in Figjsjf)), giving 7 = 12/14 = 0.85 the core appears 
(shown in red in Fig[8]F)), and the basin boundaries become fractal, see Fig[8j3)-d). 

H. Polynomial continuous-time complexity for a fixed number of formulae. We now show 
that Eq (2) of the main text implies the stronger result of having polynomial continuous - time 
complexity even for leaving fixed number (not fraction!) of formulae unsolved in the limit 

^ oc. For a given N and M = aN there are a total of 

/c-SAT formulae. If c is a small integer constant (independent of N), setting p = c/ oi*^^ {N) in (2) 
of the main text yields 

t{p,N) = b-^N^ [ln(r/c) + \nQ^^\N) 
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Figure 8: 3-XORSAT. a) basins of attraction to solutions for a formula given by the hypergraph in 
e), which has no hyperloops (core), b) basins of attraction to solutions for the formula in f), which 
has a core (shown in red), c) is the magnification of the small black rectangle from b), and d) is 
a magnification of a small black rectangle fron^^). In the hypergraphs e) and f) the "triangles" 
represent hyperedges/parity checks (k = 3) and the signs correspond to yj-s. 
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However, for fixed a. In 9^^ (TV) = a(A: - 1)7V In TV + 0(ln TV) as TV ^ oo, which implies that 

t{p,N) -7V^+i In TV, 

i.e., still showing polynomial CT complexity. Thus, assuming that (2) holds for all A^, the fact that 
only a constant number of problems may not be solved in polynomial time for ^ oo indicates 
that the algorithm runs in polynomial continuous-time on almost all hard fc-SAT instances (the 
probability to find an instance not solvable in polynomial time by this solver is zero in the N ^ oo 
limit). 



I. Computational complexity in the number of discrete steps. Since the numerical integration 
happens on a Turing machine, the true continuous trajectory is being approximated by the Runge- 
Kutta (RK) algorithm |[38l with sufficiently many discrete points, lying close to it, within a tube 
of preset diameter e (see Fig. [9]). When we monitor the fraction of formulae left unsolved p as 
function of the number of discretization steps Ugtep taken by the algorithm for hard SAT formu- 
lae (from the frozen phase), we find that p{nstep) has a power-law decay, well approximated by 



pi^step) = u{v + nstep)~^, scc Fig. 10 i. The exponent 77 also has a power-law A^-dependence 
(see Fig. IO 3): r]{N) = dN~^ with 6 2^ 1.09 ^ 1 (the A^-dependence of u and v are weak). This 
implies an exponential behavior for the number of time steps Ugtepip^ N) needed to miss solving 
only a ^-th fraction of the formulae: 

nstep{p,N) = e d p (24) 



showing exponential time-complexity for the discretized algorithm ran by a digital computer (Tur- 
ing machine). For easy formulae, however, (such as those drawn at random for a = 3), p{nstep) 
has an exponential decay, just as p{t), implying polynomial complexity for the discretized algo- 
rithm as well, see the inset of Fig[TO^. 




Figure 9: Approximating the continuous-time trajectory. The RK algorithm must compute 
a large number of discrete points in high curvature regions in order to stay within a prescribed 
precision e. 
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Figure 10: Discrete time complexity, a: Fraction of problems p{nstep) left unsolved after Ustep 
discretization steps, at o; = 4.25, A: = 3 for different system sizes N = 20, 30, 40, 50, 60, 80, 100, 
125, 150 (different colors). Averages were done over 10^ instances for each N, except for N = 
150 where 3 X 10^ instances were used. The black continuous lines are ^(nstep) = u{v + nstep)~^ - 
The inset (log-lin plot) shows p{nstep) from 10^ instances with N = 5000 at a = 3. In this easy- 
SAT region p{nstep) shows a similar exponential decay as p{t) shown in Fig3a, main text, b: The 
exponent 77 follows: r]{N) = dN~^, 5 ^ 1.09. 

J. Dependence on discretization error. The continuous-time variable t and the decay of p as 
function of t is only weakly dependent on the discretization error e by the RK solver, with a slight 
shift towards even cleaner exponentials when lowering e, see Fig. [TT^ . 

As stated in the main text, the polynomial CT complexity is not the consequence of a log- 
transformation on the time variable t. The scaling of the length L of the continuous-time tra- 
jectory measured from the initial point until it finds a solution also scales polynomially with A^. 
However, this measure is more sensitively affected by the discretization error e. This is because, 
as illustrated in Fig. |9) the length computed as the sum of the lengths of the straight segments 
between discretization steps overestimates the analog trajectory length (continuous red line). Fig. 
[TTJ3 shows the fraction of formulae p left unsolved by trajectories of length not longer than L 
(measured in T-Ln), as function of L for different discretization errors e. The convergence to clean 
exponentials is evident as e is lowered. 

K. Analogy with fluid turbulence Continuous-time processes are common in nature, from fluid 
flows to information processing in the brain; even our perception of time is arguably of analog 
nature. While transient chaos can appear in many dynamical systems, long chaotic transients 
typically occur in a parameter region preceding the permanently chaotic regime. If the average 
lifetime of transients in a dynamical system depends on an extensive parameter of the system, 
supertransients (or superpersistent chaotic transients) may appear |[39l|40|. The only experiments 
confirming the existence of such supertransients are those done in long-pipe flow experiments. 
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100 f 200 300 




L 

Figure 11: Sensitivity to discretization error, a: Fraction of problems p{t) left unsolved by 
continuous-time t, at a = 4.25, A: = 3 for system sizes N = 20, 30, 40, 50 (different colors) and 
for different error parameters of the adaptive Runge-Kutta method ll38]| : e = 0.0001, 0.001, 0.01. 
Here e is the maximal relative error allowed during the RK integration. Statistics was done on 
2.5 X 10^ instances for each N. b: The fraction of problems p{L) left unsolved by trajectories of 
length at most L (measured in Hn), for the same problems as in a. 
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found solution: t=8.2 
"2 4 6 




found solution: t=96.4 — " 
20 40 60 80 100 



20 40 60 80 100 





50 100 150 200 250 




50 100 150 200 250 



Figure 12: Time evolution of variables and energy functions. Time evolution of three variables 
Si{t) (different colors) for a formula with = 3, = 200 and a: o; = 3, characterized by 
rapid and straightforward convergence; in b: at o: = 4.25, presenting much longer, and chaotic 
trajectories, and c: for an unsatisfiable formula at a = 4.5. d: Time series of six different auxiliary 
variables a^(t) shown with different colors for the same problem as in b. Time series of E{s) 
and V^(s, a) as function of the continuous time t, for two different problems (red and black) with 
= 3, = 200 ate: a = 4.25 and f: for two unsatisfiable formulae at a = 4.5. 
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50 100 le+3 le+4 le+5 

f (analog) n^tep (discrete) 



Figure 13: Exponential decay for a single instance. For single 3-SAT instance with N = 40, 
at a = 4.25 we start the dynamics from 60000 different random initial conditions, a: Similarly to 
Fig 3a of the main text the fraction of trajectories q{t), which did not find the solution by analog- 
time t shows an exponential decay as function of t. b: The fraction of trajectories, q{nstep) still 
searching for the solution after rigtep discretization steps, however, shows a power-law behavior. 
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Figure 14: Intermittent behaviour, a: The velocity \ds/dt\ and the acceleration \(Ps/dt'^\ along 
the trajectories in the hard-SAT phase (N = 100, a = 4.25, 3-SAT) show intermittency as func- 
tion of time similar to turbulent flows [43] . b: Same as (a) on log-linear scale. 
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Here, turbulent/chaotic behavior appears before the flow enters its laminar phase (paraboHc veloc- 
ity profile), which is its only asymptotic attractor ||4T1|42|. One can think of fluid turbulence in this 
case as nature's search for equilibrium, practically solving a hard global optimization problem. Fig 
[14] perhaps takes this analogy further: it shows the fluctuations of the instantaneous velocity and 
acceleration for a typical trajectory for a hard formula from the frozen region as function of time, 
revealing intermittent behavior, typically found in turbulence, see Fig 1 of the paper by Meneveau 
and Sreenivasan 143J. 
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